!-----------------------------------------------------------------------
! $Id: fill_holes.F90 8738 2018-07-19 19:58:53Z bmg2@uwm.edu $
!===============================================================================
module fill_holes

  implicit none

  public :: fill_holes_driver, &
            fill_holes_vertical, &
            hole_filling_hm_one_lev, &
            fill_holes_hydromet, &
            fill_holes_wv, &
            vertical_avg, &
            vertical_integral, &
            clip_hydromet_conc_mvr, &
            setup_stats_indices

  private :: fill_holes_multiplicative

  private ! Set Default Scope

  contains

  !=============================================================================
  subroutine fill_holes_vertical( num_draw_pts, threshold, field_grid, &
                                  rho_ds, rho_ds_zm, &
                                  field )

    ! Description:
    ! This subroutine clips values of 'field' that are below 'threshold' as much
    ! as possible (i.e. "fills holes"), but conserves the total integrated mass
    ! of 'field'.  This prevents clipping from acting as a spurious source.
    !
    ! Mass is conserved by reducing the clipped field everywhere by a constant
    ! multiplicative coefficient.
    !
    ! This subroutine does not guarantee that the clipped field will exceed
    ! threshold everywhere; blunt clipping is needed for that.

    ! References:
    !   ``Numerical Methods for Wave Equations in Geophysical Fluid
    !     Dynamics'', Durran (1999), p. 292.
    !-----------------------------------------------------------------------

    use grid_class, only: & 
       gr ! Variable

    use clubb_precision, only: &
      core_rknd ! Variable(s)

    implicit none

    ! Input variables
    integer, intent(in) :: & 
      num_draw_pts  ! The number of points on either side of the hole;
               ! Mass is drawn from these points to fill the hole.  []

    real( kind = core_rknd ), intent(in) :: & 
      threshold  ! A threshold (e.g. w_tol*w_tol) below which field must not
                 ! fall                           [Units vary; same as field]

    character(len=2), intent(in) :: & 
      field_grid ! The grid of the field, either zt or zm

    real( kind = core_rknd ), dimension(gr%nz), intent(in) ::  & 
      rho_ds,    & ! Dry, static density on thermodynamic levels    [kg/m^3]
      rho_ds_zm    ! Dry, static density on momentum levels         [kg/m^3]

    ! Input/Output variable
    real( kind = core_rknd ), dimension(gr%nz), intent(inout) :: & 
      field  ! The field (e.g. wp2) that contains holes [Units same as threshold]

    ! Local Variables
    integer :: & 
      k,             & ! Loop index for absolute grid level              []
      begin_idx,     & ! Lower grid level of local hole-filling range    []
      end_idx,       & ! Upper grid level of local hole-filling range    []
      upper_hf_level   ! Upper grid level of global hole-filling range   []

    !-----------------------------------------------------------------------

    ! Check whether any holes exist in the entire profile.
    ! The lowest level (k=1) should not be included, as the hole-filling scheme
    ! should not alter the set value of 'field' at the surface (for momentum
    ! level variables), or consider the value of 'field' at a level below the
    ! surface (for thermodynamic level variables).  For momentum level variables
    ! only, the hole-filling scheme should not alter the set value of 'field' at
    ! the upper boundary level (k=gr%nz).

    if ( field_grid == "zt" ) then
      ! 'field' is on the zt (thermodynamic level) grid
      upper_hf_level = gr%nz
    elseif ( field_grid == "zm" )  then
      ! 'field' is on the zm (momentum level) grid
      upper_hf_level = gr%nz-1
    endif

    if ( any( field( 2:upper_hf_level ) < threshold ) ) then

      ! Make one pass up the profile, filling holes as much as we can using
      ! nearby mass.
      ! The lowest level (k=1) should not be included in the loop, as the
      ! hole-filling scheme should not alter the set value of 'field' at the
      ! surface (for momentum level variables), or consider the value of
      ! 'field' at a level below the surface (for thermodynamic level
      ! variables).  For momentum level variables only, the hole-filling scheme
      ! should not alter the set value of 'field' at the upper boundary
      ! level (k=gr%nz).
      do k = 2+num_draw_pts, upper_hf_level-num_draw_pts, 1

        begin_idx = k - num_draw_pts
        end_idx   = k + num_draw_pts

        if ( any( field( begin_idx:end_idx ) < threshold ) ) then

          ! 'field' is on the zt (thermodynamic level) grid
          if ( field_grid == "zt" ) then
            call fill_holes_multiplicative &
                    ( begin_idx, end_idx, threshold, &
                      rho_ds(begin_idx:end_idx), gr%dzt(begin_idx:end_idx), &
                      field(begin_idx:end_idx) )
                      
          ! 'field' is on the zm (momentum level) grid
          elseif ( field_grid == "zm" )  then
            call fill_holes_multiplicative &
                    ( begin_idx, end_idx, threshold, &
                      rho_ds_zm(begin_idx:end_idx), gr%dzm(begin_idx:end_idx), &
                      field(begin_idx:end_idx) )
          endif

        endif

      enddo

      ! Fill holes globally, to maximize the chance that all holes are filled.
      ! The lowest level (k=1) should not be included, as the hole-filling
      ! scheme should not alter the set value of 'field' at the surface (for
      ! momentum level variables), or consider the value of 'field' at a level
      ! below the surface (for thermodynamic level variables).  For momentum
      ! level variables only, the hole-filling scheme should not alter the set
      ! value of 'field' at the upper boundary level (k=gr%nz).
      if ( any( field( 2:upper_hf_level ) < threshold ) ) then

        ! 'field' is on the zt (thermodynamic level) grid
        if ( field_grid == "zt" ) then
          call fill_holes_multiplicative &
                 ( 2, upper_hf_level, threshold, &
                   rho_ds(2:upper_hf_level), gr%dzt(2:upper_hf_level), &
                   field(2:upper_hf_level) )
                   
        ! 'field' is on the zm (momentum level) grid
        elseif ( field_grid == "zm" )  then
            call fill_holes_multiplicative &
                 ( 2, upper_hf_level, threshold, &
                   rho_ds_zm(2:upper_hf_level), gr%dzm(2:upper_hf_level), &
                   field(2:upper_hf_level) )
        endif

      endif

    endif  ! End overall check for existence of holes

    return

  end subroutine fill_holes_vertical

  !=============================================================================
  subroutine fill_holes_multiplicative &
                 ( begin_idx, end_idx, threshold, &
                   rho, dz, &
                   field )

    ! Description:
    ! This subroutine clips values of 'field' that are below 'threshold' as much
    ! as possible (i.e. "fills holes"), but conserves the total integrated mass
    ! of 'field'.  This prevents clipping from acting as a spurious source.
    !
    ! Mass is conserved by reducing the clipped field everywhere by a constant
    ! multiplicative coefficient.
    !
    ! This subroutine does not guarantee that the clipped field will exceed
    ! threshold everywhere; blunt clipping is needed for that.

    ! References:
    ! ``Numerical Methods for Wave Equations in Geophysical Fluid
    ! Dynamics", Durran (1999), p. 292.
    !-----------------------------------------------------------------------

    use constants_clubb, only: &
      eps

    use clubb_precision, only: &
      core_rknd ! Variable(s)

    implicit none

    ! Input variables
    integer, intent(in) :: & 
      begin_idx, & ! The beginning index (e.g. k=2) of the range of hole-filling 
      end_idx      ! The end index (e.g. k=gr%nz) of the range of hole-filling

    real( kind = core_rknd ), intent(in) :: & 
      threshold  ! A threshold (e.g. w_tol*w_tol) below which field must not fall
                 !                              [Units vary; same as field]

    real( kind = core_rknd ), dimension(end_idx-begin_idx+1), intent(in) ::  & 
      rho,     &  ! Dry, static density on either thermodynamic or momentum levels   [kg/m^3]
      dz          ! Reciprocal of thermodynamic or momentum level thickness depending on whether
                  ! we're on zt or zm grid.

    ! Input/Output variable
    real( kind = core_rknd ), dimension(end_idx-begin_idx+1), intent(inout) ::  & 
      field  ! The field (e.g. wp2) that contains holes
             !                                  [Units same as threshold]

    ! Local Variables
    real( kind = core_rknd ), dimension(end_idx-begin_idx+1)  ::  & 
      field_clipped  ! The raw field (e.g. wp2) that contains no holes
                     !                          [Units same as threshold]

    real( kind = core_rknd ) ::  & 
      field_avg,  &        ! Vertical average of field [Units of field]
      field_clipped_avg, & ! Vertical average of clipped field [Units of field]
      mass_fraction        ! Coefficient that multiplies clipped field
                           ! in order to conserve mass.                      []

    !-----------------------------------------------------------------------

    ! Compute the field's vertical average, which we must conserve.
    field_avg = vertical_avg( (end_idx-begin_idx+1), rho, &
                                  field, dz )

    ! Clip small or negative values from field.
    if ( field_avg >= threshold ) then
      ! We know we can fill in holes completely
      field_clipped = max( threshold, field )
    else
      ! We can only fill in holes partly;
      ! to do so, we remove all mass above threshold.
      field_clipped = min( threshold, field )
    endif

    ! Compute the clipped field's vertical integral.
    ! clipped_total_mass >= original_total_mass
    field_clipped_avg = vertical_avg( (end_idx-begin_idx+1), rho, &
                                      field_clipped, dz )

    ! If the difference between the field_clipped_avg and the threshold is so
    ! small that it falls within numerical round-off, return to the parent
    ! subroutine without altering the field in order to avoid divide-by-zero
    ! error.
    if ( abs(field_clipped_avg-threshold) <= abs(field_clipped_avg+threshold)*eps/2) then
      return
    endif

    ! Compute coefficient that makes the clipped field have the same mass as the
    ! original field.  We should always have mass_fraction > 0.
    mass_fraction = ( field_avg - threshold ) / & 
                          ( field_clipped_avg - threshold )

    ! Output normalized, filled field
    field = mass_fraction * ( field_clipped - threshold )  & 
                 + threshold

    return

  end subroutine fill_holes_multiplicative

  !=============================================================================
  function vertical_avg( total_idx, rho_ds, field, dz )

    ! Description:
    ! Computes the density-weighted vertical average of a field.
    !
    ! The average value of a function, f, over a set domain, [a,b], is
    ! calculated by the equation:
    !
    ! f_avg = ( INT(a:b) f*g ) / ( INT(a:b) g );
    !
    ! as long as f is continous and g is nonnegative and integrable.  Therefore,
    ! the density-weighted (by dry, static, base-static density) vertical
    ! average value of any model field, x, is calculated by the equation:
    !
    ! x_avg|_z = ( INT(z_bot:z_top) x rho_ds dz )
    !            / ( INT(z_bot:z_top) rho_ds dz );
    !
    ! where z_bot is the bottom of the vertical domain, and z_top is the top of
    ! the vertical domain.
    !
    ! This calculation is done slightly differently depending on whether x is a
    ! thermodynamic-level or a momentum-level variable.
    !
    ! Thermodynamic-level computation:
    
    !
    ! For numerical purposes, INT(z_bot:z_top) x rho_ds dz, which is the
    ! numerator integral, is calculated as:
    !
    ! SUM(k_bot:k_top) x(k) rho_ds(k) delta_z(k);
    !
    ! where k is the index of the given thermodynamic level, x and rho_ds are
    ! both thermodynamic-level variables, and delta_z(k) = zm(k) - zm(k-1).  The
    ! indices k_bot and k_top are the indices of the respective lower and upper
    ! thermodynamic levels involved in the integration.
    !
    ! Likewise, INT(z_bot:z_top) rho_ds dz, which is the denominator integral,
    ! is calculated as:
    !
    ! SUM(k_bot:k_top) rho_ds(k) delta_z(k).
    !
    ! The first (k=1) thermodynamic level is below ground (or below the
    ! official lower boundary at the first momentum level), so it should not
    ! count in a vertical average, whether that vertical average is used for
    ! the hole-filling scheme or for statistical purposes. Begin no lower
    ! than level k=2, which is the first thermodynamic level above ground (or
    ! above the model lower boundary).
    !
    ! For cases where hole-filling over the entire (global) vertical domain
    ! is desired, or where statistics over the entire (global) vertical
    ! domain are desired, the lower (thermodynamic-level) index of k = 2 and
    ! the upper (thermodynamic-level) index of k = gr%nz, means that the
    ! overall vertical domain will be gr%zm(gr%nz) - gr%zm(1).
    !
    !
    ! Momentum-level computation:
    !
    ! For numerical purposes, INT(z_bot:z_top) x rho_ds dz, which is the
    ! numerator integral, is calculated as:
    !
    ! SUM(k_bot:k_top) x(k) rho_ds(k) delta_z(k);
    !
    ! where k is the index of the given momentum level, x and rho_ds are both
    ! momentum-level variables, and delta_z(k) = zt(k+1) - zt(k).  The indices
    ! k_bot and k_top are the indices of the respective lower and upper momentum
    ! levels involved in the integration.
    !
    ! Likewise, INT(z_bot:z_top) rho_ds dz, which is the denominator integral,
    ! is calculated as:
    !
    ! SUM(k_bot:k_top) rho_ds(k) delta_z(k).
    !
    ! The first (k=1) momentum level is right at ground level (or right at
    ! the official lower boundary).  The momentum level variables that call
    ! the hole-filling scheme have set values at the surface (or lower
    ! boundary), and those set values should not be changed.  Therefore, the
    ! vertical average (for purposes of hole-filling) should not include the
    ! surface level (or lower boundary level).  For hole-filling purposes,
    ! begin no lower than level k=2, which is the second momentum level above
    ! ground (or above the model lower boundary).  Likewise, the value at the
    ! model upper boundary (k=gr%nz) is also set for momentum level
    ! variables.  That value should also not be changed.
    !
    ! However, this function is also used to keep track (for statistical
    ! purposes) of the vertical average of certain variables.  In that case,
    ! the vertical average needs to be taken over the entire vertical domain
    ! (level 1 to level gr%nz).
    !
    !
    ! In both the thermodynamic-level computation and the momentum-level
    ! computation, the numerator integral is divided by the denominator integral
    ! in order to find the average value (over the vertical domain) of x.

    ! References:
    ! None
    !-----------------------------------------------------------------------

    use clubb_precision, only: &
      core_rknd ! Variable(s)

    implicit none

    ! Input variables
    integer, intent(in) :: & 
      total_idx ! The total numer of indices within the range of averaging

    real( kind = core_rknd ), dimension(total_idx), intent(in) ::  &
      rho_ds, & ! Dry, static density on either thermodynamic or momentum levels    [kg/m^3]
      field,  & ! The field (e.g. wp2) to be vertically averaged                    [Units vary]
      dz  ! Reciprocal of thermodynamic or momentum level thickness           [1/m]
                ! depending on whether we're on zt or zm grid.
    ! Note:  The rho_ds and field points need to be arranged from
    !        lowest to highest in altitude, with rho_ds(1) and
    !        field(1) actually their respective values at level k = 1.

    ! Output variable
    real( kind = core_rknd ) :: & 
      vertical_avg  ! Vertical average of field    [Units of field]

    ! Local variables
    real( kind = core_rknd ) :: & 
      numer_integral, & ! Integral in the numerator (see description)
      denom_integral    ! Integral in the denominator (see description)
      

    integer :: k

    !-----------------------------------------------------------------------
    
    ! Initialize variable
    numer_integral = 0.0_core_rknd
    denom_integral = 0.0_core_rknd

    ! Compute the numerator and denominator integral.
    ! Multiply rho_ds at level k by the level thickness
    ! at level k.  Then, sum over all vertical levels.
    do k=1, total_idx

        numer_integral = numer_integral + rho_ds(k) * dz(k) * field(k)
        denom_integral = denom_integral + rho_ds(k) * dz(k)

    end do

    ! Find the vertical average of 'field'.
    vertical_avg = numer_integral / denom_integral

    return
  end function vertical_avg

  !=============================================================================
  pure function vertical_integral( total_idx, rho_ds, &
                                       field, dz )

    ! Description:
    ! Computes the vertical integral. rho_ds, field, and dz must all be
    ! of size total_idx and should all start at the same index.
    ! 
    
    ! References:
    ! None
    !-----------------------------------------------------------------------

    use clubb_precision, only: &
      core_rknd ! Variable(s)

    implicit none

    ! Input variables
    integer, intent(in) :: & 
      total_idx  ! The total numer of indices within the range of averaging

    real( kind = core_rknd ), dimension(total_idx), intent(in) ::  &
      rho_ds,  & ! Dry, static density                   [kg/m^3]
      field,   & ! The field to be vertically averaged   [Units vary]
      dz         ! Level thickness                       [1/m]
    ! Note:  The rho_ds and field points need to be arranged from
    !        lowest to highest in altitude, with rho_ds(1) and
    !        field(1) actually their respective values at level k = begin_idx.

    ! Local variables
    real( kind = core_rknd ) :: &
      vertical_integral ! Integral in the numerator (see description)

    !-----------------------------------------------------------------------

    !  Assertion checks: that begin_idx <= gr%nz - 1
    !                    that end_idx   >= 2
    !                    that begin_idx <= end_idx


    ! Initializing vertical_integral to avoid a compiler warning.
    vertical_integral = 0.0_core_rknd

    ! Compute the integral.
    ! Multiply the field at level k by rho_ds at level k and by
    ! the level thickness at level k.  Then, sum over all vertical levels.
    ! Note:  The values of the field and rho_ds are passed into this function
    !        so that field(1) and rho_ds(1) are actually the field and rho_ds
    !        at level k_start.
    vertical_integral = sum( field * rho_ds * dz )

    !print *, vertical_integral

    return
  end function vertical_integral

!===============================================================================

  subroutine hole_filling_hm_one_lev( num_hm_fill, hm_one_lev, & ! Intent(in)
                                   hm_one_lev_filled ) ! Intent(out)

  ! Description:
  ! Fills holes between same-phase (i.e. either liquid or frozen) hydrometeors for
  ! one height level.
  !
  ! Warning: Do not input hydrometeors of different phases, e.g. liquid and frozen.
  ! Otherwise heat will not be conserved.
  !
  ! References:
  !
  ! None
  !-----------------------------------------------------------------------

    use constants_clubb, only: &
        one, & ! Variable(s)
        zero, &
        eps

    use clubb_precision, only: &
        core_rknd ! Variable(s)

    use error_code, only: &
        clubb_at_least_debug_level  ! Procedure

    implicit none

    ! Input Variables
    integer, intent(in) :: num_hm_fill ! number of hydrometeors involved

    real(kind = core_rknd), dimension(num_hm_fill), intent(in) :: hm_one_lev

    ! Output Variables
    real(kind = core_rknd), dimension(num_hm_fill), intent(out) :: hm_one_lev_filled

    ! Local Variables
    integer :: num_neg_hm ! number of holes

    real(kind = core_rknd) :: &
      total_hole, & ! Size of the hole ( missing mass, less than 0 )
      total_mass    ! Total mass to fill the hole
      ! total mass of water substance = total_mass + total_hole

    integer :: i ! loop iterator

  !-----------------------------------------------------------------------

    !----- Begin Code -----

    ! Initialization
    hm_one_lev_filled = 0._core_rknd
    total_hole = 0._core_rknd
    total_mass = 0._core_rknd
    num_neg_hm = 0

    ! Determine the total size of the hole and the number of neg. hydrometeors
    ! and the total mass of hole filling material
    do i=1, num_hm_fill
!       print *, "hm_one_lev(",i,") = ", hm_one_lev(i)
       if ( hm_one_lev(i) < zero ) then
          total_hole = total_hole + hm_one_lev(i) ! less than zero
          num_neg_hm = num_neg_hm + 1
       else
          total_mass = total_mass + hm_one_lev(i)
       endif

    enddo

!    print *, "total_hole = ", total_hole
!    print *, "total_mass = ", total_mass
!    print *, "num_neg_hm = ", num_neg_hm

    ! There is no water substance at all to fill the hole
    if ( abs(total_mass) < eps ) then

       if ( clubb_at_least_debug_level( 2 ) ) then
          print *, "Warning: One-level hole filling was not successful! total_mass ~= 0"
       endif

       hm_one_lev_filled = hm_one_lev

       return
    endif

    ! Fill the holes and adjust the remaining quantities:
    ! hm_filled(i) = 0, if hm(i) < 0
    ! or
    ! hm_filled(i) = (1 + total_hole/total_mass)*hm(i), if hm(i) > 0
    do i=1, num_hm_fill

       ! if there is not enough material, fill the holes partially with all the material available
       if ( abs(total_hole) > total_mass ) then

          if ( clubb_at_least_debug_level( 2 ) ) then
             print *, "Warning: One-level hole filling was not able to fill holes completely!" // &
                      " The holes were filled partially. |total_hole| > total_mass"
          endif

          hm_one_lev_filled(i) = min(hm_one_lev(i), zero) * ( one + total_mass / total_hole )

       else ! fill holes completely
          hm_one_lev_filled(i) = max(hm_one_lev(i), zero) * ( one + total_hole / total_mass )

       endif

    enddo

    ! Assertion checks (water substance conservation, non-negativity)
    if ( clubb_at_least_debug_level( 2 ) ) then

       if ( abs(sum( hm_one_lev ) - sum(hm_one_lev_filled)) > &
            abs(sum( hm_one_lev ) + sum(hm_one_lev_filled)) * eps/2 ) then
          print *, "Warning: Hole filling was not conservative!"
       endif

       if ( any( hm_one_lev_filled < zero ) ) then
          print *, "Warning: Hole filling failed! A hole could not be filled."
       endif

    endif

    return

  end subroutine hole_filling_hm_one_lev
  !-----------------------------------------------------------------------

  !-----------------------------------------------------------------------
  subroutine fill_holes_hydromet( nz, hydromet_dim, hydromet, & ! Intent(in)
                                  hydromet_filled ) ! Intent(out)

  ! Description:
  ! Fills holes between same-phase hydrometeors(i.e. for frozen hydrometeors).
  ! The hole filling conserves water substance between all same-phase (frozen or liquid)
  ! hydrometeors at each height level.
  !
  ! Attention: The hole filling for the liquid phase hydrometeors is not yet implemented
  !
  ! Attention: l_frozen_hm and l_mix_rat_hm need to be set up before this subroutine is called!
  !
  ! References:
  !
  ! None
  !-----------------------------------------------------------------------

    use clubb_precision, only: &
        core_rknd

    use array_index, only: &
        l_frozen_hm, & ! Variable(s)
        l_mix_rat_hm

    use constants_clubb, only: &
        zero

    implicit none

    ! Input Variables
    integer, intent(in) :: hydromet_dim, nz

    real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(in) :: &
      hydromet

    ! Output Variables
    real( kind = core_rknd ), dimension(nz,hydromet_dim), intent(out) :: &
      hydromet_filled

    ! Local Variables
    integer :: i,j ! Loop iterators

    integer :: num_frozen_hm ! Number of frozen hydrometeor mixing ratios

    real( kind = core_rknd ), dimension(:,:), allocatable :: &
      hydromet_frozen,       & ! Frozen hydrometeor mixing ratios
      hydromet_frozen_filled   ! Frozen hydrometeor mixing ratios after hole filling

  !-----------------------------------------------------------------------

    !----- Begin Code -----

    ! Determine the number of frozen hydrometeor mixing ratios
    num_frozen_hm = 0
    do i=1,hydromet_dim
       if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
          num_frozen_hm = num_frozen_hm + 1
       endif
    enddo

    ! Allocation
    allocate( hydromet_frozen(nz,num_frozen_hm) )
    allocate( hydromet_frozen_filled(nz,num_frozen_hm) )

    ! Determine frozen hydrometeor mixing ratios
    j = 1
    do i = 1,hydromet_dim
       if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
          hydromet_frozen(:,j) = hydromet(:,i)
          j = j+1
       endif
    enddo

    ! Fill holes for the frozen hydrometeors
    do i=1,nz
       if ( any( hydromet_frozen(i,:) < zero ) ) then
          call hole_filling_hm_one_lev( num_frozen_hm, hydromet_frozen(i,:), & ! Intent(in)
                                     hydromet_frozen_filled(i,:) ) ! Intent(out)
       else
          hydromet_frozen_filled(i,:) = hydromet_frozen(i,:)
       endif
    enddo

    ! Setup the filled hydromet array
    j = 1
    do i=1, hydromet_dim
       if ( l_frozen_hm(i) .and. l_mix_rat_hm(i) ) then
          hydromet_filled(:,i) = hydromet_frozen_filled(:,j)
          j = j+1
       else
          hydromet_filled(:,i) = hydromet(:,i)
       endif
    enddo

    !!! Here we could do the same hole filling for all the liquid phase hydrometeors

    return
  end subroutine fill_holes_hydromet
  !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
  subroutine fill_holes_wv( nz, dt, exner, hydromet_name, & ! Intent(in)
                            rvm_mc, thlm_mc, hydromet )! Intent(inout)

  ! Description:
  ! Fills holes using the cloud water mixing ratio from the current height level.
  !
  ! References:
  !
  ! None
  !-----------------------------------------------------------------------

    use clubb_precision, only: &
        core_rknd

    use constants_clubb, only: &
        zero_threshold, &
        Lv, &
        Ls, &
        Cp

    implicit none

    ! Input Variables
    integer, intent(in) :: nz

    real( kind = core_rknd ), intent(in) ::  &
      dt           ! Timestep         [s]

    character(len=10), intent(in) :: hydromet_name

    real( kind = core_rknd ), dimension(nz), intent(in) :: &
      exner        ! Exner function                            [-]

    ! Input/Output Variables
    real( kind = core_rknd ), dimension(nz), intent(inout) :: &
      hydromet, &  ! Hydrometeor array                         [units vary]
      rvm_mc, &
      thlm_mc

    ! Local Variables
    integer :: k ! Loop iterator

    real( kind = core_rknd ) :: rvm_clip_tndcy
  !-----------------------------------------------------------------------

    !----- Begin Code -----

    do k = 2, nz, 1

       if ( hydromet(k) < zero_threshold ) then

          ! Set rvm_clip_tndcy to the time tendency applied to vapor and removed
          ! from the hydrometeor.
          rvm_clip_tndcy = hydromet(k) / dt

          ! Adjust the tendency rvm_mc accordingly
          rvm_mc(k) = rvm_mc(k) + rvm_clip_tndcy

          ! Adjust the tendency of thlm_mc according to whether the
          ! effect is an evaporation or sublimation tendency.
          select case ( trim( hydromet_name ) )
          case( "rrm" )
             thlm_mc(k) = thlm_mc(k) - rvm_clip_tndcy * ( Lv / ( Cp*exner(k) ) )
          case( "rim", "rsm", "rgm" )
             thlm_mc(k) = thlm_mc(k) - rvm_clip_tndcy * ( Ls / ( Cp*exner(k) ) )
          case default
             stop "Fatal error in microphys_driver"
          end select

          ! Set the mixing ratio to 0
          hydromet(k) = zero_threshold

       endif ! hydromet(k,i) < 0

    enddo ! k = 2..gr%nz

    return
  end subroutine fill_holes_wv
  !-----------------------------------------------------------------------

  !-----------------------------------------------------------------------
  subroutine fill_holes_driver( nz, dt, hydromet_dim,        & ! Intent(in)
                                l_fill_holes_hm,             & ! Intent(in)
                                rho_ds_zm, rho_ds_zt, exner, & ! Intent(in)
                                thlm_mc, rvm_mc, hydromet )    ! Intent(inout)

  ! Description:
  ! Fills holes between same-phase hydrometeors(i.e. for frozen hydrometeors).
  ! The hole filling conserves water substance between all same-phase (frozen or liquid)
  ! hydrometeors at each height level.
  !
  ! Attention: The hole filling for the liquid phase hydrometeors is not yet implemented
  !
  ! Attention: l_frozen_hm and l_mix_rat_hm need to be set up before this subroutine is called!
  !
  ! References:
  !
  ! None
  !-----------------------------------------------------------------------

    use grid_class, only: &
        gr  ! Variable(s)

    use clubb_precision, only: &
        core_rknd   ! Variable(s)

    use constants_clubb, only: &
        zero,            &
        zero_threshold,  &
        Lv,              &
        Ls,              &
        Cp,              &
        fstderr

    use array_index, only: &
        hydromet_list, & ! Names of the hydrometeor species
        hydromet_tol

    use array_index, only: &
        l_mix_rat_hm, & ! Variable(s)
        l_frozen_hm

    use index_mapping, only: &
        Nx2rx_hm_idx, & ! Procedure(s)
        mvr_hm_max

    use stats_type_utilities, only: &
        stat_begin_update, & ! Subroutines
        stat_end_update

    use stats_variables, only: &
        stats_zt, &  ! Variables
        l_stats_samp

    use error_code, only: &
        clubb_at_least_debug_level  ! Procedure

    implicit none

    intrinsic :: trim

    ! Input Variables
    integer, intent(in) :: hydromet_dim, nz

    logical, intent(in) :: l_fill_holes_hm

    real( kind = core_rknd ), intent(in) ::  &
      dt           ! Timestep         [s]

    real( kind = core_rknd ), dimension(nz), intent(in) :: &
      rho_ds_zm, & ! Dry, static density on momentum levels   [kg/m^3]
      rho_ds_zt    ! Dry, static density on thermo. levels    [kg/m^3]

    real( kind = core_rknd ), dimension(nz), intent(in) :: &
      exner  ! Exner function                                       [-]

    ! Input/Output Variables
    real( kind = core_rknd ), dimension(nz, hydromet_dim), intent(inout) :: &
      hydromet    ! Mean of hydrometeor fields    [units vary]

    real( kind = core_rknd ), dimension(nz), intent(inout) :: &
      rvm_mc,  & ! Microphysics contributions to vapor water           [kg/kg/s]
      thlm_mc    ! Microphysics contributions to liquid potential temp [K/s]

    ! Local Variables
    integer :: i, k ! Loop iterators

    real( kind = core_rknd ), dimension(nz, hydromet_dim) :: &
      hydromet_filled,  & ! Frozen hydrometeor mixing ratios after hole filling
      hydromet_clipped    ! Clipped mean of hydrometeor fields    [units vary]

    character( len = 10 ) :: hydromet_name

    real( kind = core_rknd ) :: &
      max_velocity    ! Maximum sedimentation velocity                     [m/s]

    integer :: ixrm_hf, ixrm_wvhf, ixrm_cl, &
               ixrm_bt, ixrm_mc

    logical :: l_hole_fill = .true.

  !-----------------------------------------------------------------------

    !----- Begin Code -----

    ! Start stats output for the _hf variables (changes in the hydromet array
    ! due to fill_holes_hydromet and fill_holes_vertical)
    if ( l_stats_samp ) then

       do i = 1, hydromet_dim

          ! Set up the stats indices for hydrometeor at index i
          call setup_stats_indices( i,                           & ! Intent(in)
                                    ixrm_bt, ixrm_hf, ixrm_wvhf, & ! Intent(inout)
                                    ixrm_cl, ixrm_mc,            & ! Intent(inout)
                                    max_velocity )                 ! Intent(inout)

          call stat_begin_update( ixrm_hf, hydromet(:,i) &
                                           / dt, stats_zt )

       enddo ! i = 1, hydromet_dim

    endif ! l_stats_samp

    ! If we're dealing with negative hydrometeors, we first try to fill the
    ! holes proportionally from other same-phase hydrometeors at each height
    ! level.
    if ( any( hydromet < zero_threshold ) .and. l_fill_holes_hm ) then

       call fill_holes_hydromet( nz, hydromet_dim, hydromet, & ! Intent(in)
                                 hydromet_filled ) ! Intent(out)

       hydromet = hydromet_filled

    endif ! any( hydromet < zero ) .and. l_fill_holes_hm

    hydromet_filled = zero

    do i = 1, hydromet_dim

      ! Set up the stats indices for hydrometeor at index i
      call setup_stats_indices( i,                           & ! Intent(in)
                                ixrm_bt, ixrm_hf, ixrm_wvhf, & ! Intent(inout)
                                ixrm_cl, ixrm_mc,            & ! Intent(inout)
                                max_velocity )                 ! Intent(inout)

      hydromet_name = hydromet_list(i)

      ! Print warning message if any hydrometeor species has a value < 0.
      if ( clubb_at_least_debug_level( 1 ) ) then
         if ( any( hydromet(:,i) < zero_threshold ) ) then

            do k = 1, nz
               if ( hydromet(k,i) < zero_threshold ) then
                  write(fstderr,*) trim( hydromet_name ) //" < ", &
                                   zero_threshold, &
                                   " in fill_holes_driver at k= ", k
               endif ! hydromet(k,i) < 0
            enddo ! k = 1, nz

         endif ! hydromet(:,i) < 0       
      endif ! clubb_at_least_debug_level( 1 )


      ! Store the previous value of the hydrometeor for the effect of the
      ! hole-filling scheme.
!      if ( l_stats_samp ) then
!         call stat_begin_update( ixrm_hf, hydromet(:,i) &
!                                          / dt, stats_zt )
!      endif

      ! If we're dealing with a mixing ratio and hole filling is enabled,
      ! then we apply the hole filling algorithm
      if ( any( hydromet(:,i) < zero_threshold ) ) then

         if ( hydromet_name(1:1) == "r" .and. l_hole_fill ) then

            ! Apply the hole filling algorithm
            call fill_holes_vertical( 2, zero_threshold, "zt", &
                                      rho_ds_zt, rho_ds_zm, &
                                      hydromet(:,i) )

         endif ! Variable is a mixing ratio and l_hole_fill is true

      endif ! hydromet(:,i) < 0

      ! Enter the new value of the hydrometeor for the effect of the
      ! hole-filling scheme.
      if ( l_stats_samp ) then
         call stat_end_update( ixrm_hf, hydromet(:,i) &
                                        / dt, stats_zt )
      endif

      ! Store the previous value of the hydrometeor for the effect of the water
      ! vapor hole-filling scheme.
      if ( l_stats_samp ) then
         call stat_begin_update( ixrm_wvhf, hydromet(:,i) &
                                            / dt, stats_zt )
      endif

      if ( any( hydromet(:,i) < zero_threshold ) ) then

         if ( hydromet_name(1:1) == "r" .and. l_hole_fill ) then

            ! If the hole filling algorithm failed, then we attempt to fill
            ! the missing mass with water vapor mixing ratio.
            ! We noticed this is needed for ASEX A209, particularly if Latin
            ! hypercube sampling is enabled.  -dschanen 11 Nov 2010
            call fill_holes_wv( nz, dt, exner, hydromet_name, & ! Intent(in)
                                rvm_mc, thlm_mc, hydromet(:,i) )   ! Intent(out)

         endif ! Variable is a mixing ratio and l_hole_fill is true

      endif ! hydromet(:,i) < 0

      ! Enter the new value of the hydrometeor for the effect of the water vapor
      ! hole-filling scheme.
      if ( l_stats_samp ) then
         call stat_end_update( ixrm_wvhf, hydromet(:,i) &
                                          / dt, stats_zt )
      endif

      ! Clipping for hydrometeor mixing ratios.
      if ( l_mix_rat_hm(i) ) then

         ! Store the previous value of the hydrometeor for the effect of
         ! clipping.
         if ( l_stats_samp ) then
            call stat_begin_update( ixrm_cl, &
                                    hydromet(:,i) &
                                    / dt, &
                                    stats_zt )
         endif

         if ( any( hydromet(:,i) < zero_threshold ) ) then

            ! Clip any remaining negative values of precipitating hydrometeor
            ! mixing ratios to 0.
            where ( hydromet(:,i) < zero_threshold )
               hydromet(:,i) = zero_threshold
            end where

         endif ! hydromet(:,i) < 0

         ! Eliminate very small values of mean precipitating hydrometeor mixing
         ! ratios by setting them to 0.
         do k = 2, gr%nz, 1

            if ( hydromet(k,i) <= hydromet_tol(i) ) then

               rvm_mc(k) &
               = rvm_mc(k) &
                 + ( hydromet(k,i) / dt )

               if ( .not. l_frozen_hm(i) ) then

                  ! Rain water mixing ratio
   
                  thlm_mc(k) &
                  = thlm_mc(k) &
                    - ( Lv / ( Cp * exner(k) ) ) &
                      * ( hydromet(k,i) / dt )

               else ! Frozen hydrometeor mixing ratio

                  thlm_mc(k) &
                  = thlm_mc(k) &
                    - ( Ls / ( Cp * exner(k) ) ) &
                      * ( hydromet(k,i) / dt )

               endif ! l_frozen_hm(i)

               hydromet(k,i) = zero

            endif ! hydromet(k,i) <= hydromet_tol(i)

         enddo ! k = 2, gr%nz, 1


         ! Enter the new value of the hydrometeor for the effect of clipping.
         if ( l_stats_samp ) then
            call stat_end_update( ixrm_cl, hydromet(:,i) &
                                           / dt, stats_zt )
         endif

      endif ! l_mix_rat_hm(i)

    enddo ! i = 1, hydromet_dim, 1

    ! Calculate clipping for hydrometeor concentrations.
    call clip_hydromet_conc_mvr( hydromet_dim, hydromet, & ! Intent(in)
                                 hydromet_clipped )        ! Intent(out)

    ! Clip hydrometeor concentrations and output stats.
    do i = 1, hydromet_dim

       if ( .not. l_mix_rat_hm(i) ) then

          if ( l_stats_samp ) then

             ! Set up the stats indices for hydrometeor at index i
             call setup_stats_indices( i,                           & ! In
                                       ixrm_bt, ixrm_hf, ixrm_wvhf, & ! In/Out
                                       ixrm_cl, ixrm_mc,            & ! In/Out
                                       max_velocity )                 ! In/Out

             ! Store the previous value of the hydrometeor for the effect of
             ! clipping.
             call stat_begin_update( ixrm_cl, hydromet(:,i) / dt, stats_zt )

          endif ! l_stats_samp

          ! Apply clipping of hydrometeor concentrations.
          hydromet(:,i) = hydromet_clipped(:,i)

          ! Enter the new value of the hydrometeor for the effect of clipping.
          if ( l_stats_samp ) then
             call stat_end_update( ixrm_cl, hydromet(:,i) / dt, stats_zt )
          endif

       endif ! .not. l_mix_rat_hm(i)

    enddo ! i = 1, hydromet_dim, 1


    return

  end subroutine fill_holes_driver

  !=============================================================================
  subroutine clip_hydromet_conc_mvr( hydromet_dim, hydromet, & ! Intent(in)
                                     hydromet_clipped )        ! Intent(out)

    ! Description:
    ! Increases the value of a hydrometeor concentration when it is too small,
    ! according to the hydrometeor mixing ratio (which remains unchanged) and
    ! the maximum acceptable drop or particle mean volume radius for that
    ! hydrometeor species.

    ! References:
    !-----------------------------------------------------------------------

    use grid_class, only: &
        gr    ! Variable(s)

    use constants_clubb, only: &
        pi,          & ! Variable(s)
        four_thirds, &
        one,         &
        zero,        &
        rho_lw,      &
        rho_ice

    use array_index, only: &
        l_mix_rat_hm, & ! Variable(s)
        l_frozen_hm

    use index_mapping, only: &
        Nx2rx_hm_idx, & ! Procedure(s)
        mvr_hm_max

    use clubb_precision, only: &
        core_rknd    ! Variable(s)

    implicit none

    ! Input Variables
    integer, intent(in) :: &
      hydromet_dim    ! Number of hydrometeor fields

    real( kind = core_rknd ), dimension(gr%nz,hydromet_dim), intent(in) :: &
      hydromet    ! Mean of hydrometeor fields    [units vary]

    ! Output Variable
    real( kind = core_rknd ), dimension(gr%nz,hydromet_dim), intent(out) :: &
      hydromet_clipped    ! Clipped mean of hydrometeor fields    [units vary]

    ! Local Variables
    real( kind = core_rknd ) :: &
      Nxm_min_coef    ! Coefficient for min mean value of a concentration [1/kg]

    integer :: &
      k,   & ! Vertical grid index
      idx    ! Hydrometeor species index


    ! Clipping for hydrometeor concentrations.
    do idx = 1, hydromet_dim

       if ( .not. l_mix_rat_hm(idx) ) then

          ! The field is a hydrometeor concentration.

          if ( .not. l_frozen_hm(idx) ) then

             ! Clipping for mean rain drop concentration, <Nr>.
             !
             ! When mean rain water mixing ratio, <rr>, is found at a grid
             ! level, mean rain drop concentration must be at least a minimum
             ! value so that average rain drop mean volume radius stays within
             ! an upper bound.  Otherwise, mean rain drop concentration is 0.
             ! This can also be applied to values at a sample point, rather than
             ! a grid-box mean.

             ! The minimum mean rain drop concentration is given by:
             !
             ! <Nr> = <rr> / ( (4/3) * pi * rho_lw * mvr_rain_max^3 );
             !
             ! where mvr_rain_max is the maximum acceptable average rain drop
             ! mean volume radius.

             Nxm_min_coef &
             = one / ( four_thirds * pi * rho_lw * mvr_hm_max(idx)**3 )

          else ! l_frozen_hm(idx)

             ! Clipping for mean frozen hydrometeor concentration, <Nx>, where x
             ! stands for any frozen hydrometeor species.
             !
             ! When mean frozen hydrometeor mixing ratio, <rx>, is found at a
             ! grid level, mean frozen hydrometeor concentration must be at
             ! least a minimum value so that average frozen hydrometeor mean
             ! volume radius stays within an upper bound.  Otherwise, mean
             ! frozen hydrometeor concentration is 0.  This can also be applied
             ! to values at a sample point, rather than a grid-box mean.

             ! The minimum mean frozen hydrometeor concentration is given by:
             !
             ! <Nx> = <rx> / ( (4/3) * pi * rho_ice * mvr_x_max^3 );
             !
             ! where mvr_x_max is the maximum acceptable average frozen
             ! hydrometeor mean volume radius for frozen hydrometeor species, x.

             Nxm_min_coef &
             = one / ( four_thirds * pi * rho_ice * mvr_hm_max(idx)**3 )

          endif ! .not. l_frozen_hm(idx)

          ! Loop over vertical levels and increase hydrometeor concentrations
          ! when necessary.
          do k = 2, gr%nz, 1

             if ( hydromet(k,Nx2rx_hm_idx(idx)) > zero ) then

                ! Hydrometeor mixing ratio, <rx>, is found at the grid level.
                hydromet_clipped(k,idx) &
                = max( hydromet(k,idx), &
                       Nxm_min_coef * hydromet(k,Nx2rx_hm_idx(idx)) )

             else ! <rx> = 0

                hydromet_clipped(k,idx) = zero

             endif ! hydromet(k,Nx2rx_hm_idx(idx)) > 0

          enddo ! k = 2, gr%nz, 1

          ! The lowest thermodynamic level is below the model's lower boundary.
          hydromet_clipped(1,idx) = hydromet(1,idx)

       else ! l_mix_rat_hm(idx)

          ! The field is a hydrometeor mixing ratio.
          hydromet_clipped(:,idx) = hydromet(:,idx)

       endif ! .not. l_mix_rat_hm(idx)

    enddo ! idx = 1, hydromet_dim, 1


    return

  end subroutine clip_hydromet_conc_mvr

  !-----------------------------------------------------------------------
  subroutine setup_stats_indices( ihm,                         & ! Intent(in)
                                  ixrm_bt, ixrm_hf, ixrm_wvhf, & ! Intent(inout)
                                  ixrm_cl, ixrm_mc,            & ! Intent(inout)
                                  max_velocity )                 ! Intent(inout)

  ! Description:
  !
  ! Determines the stats output indices depending on the hydrometeor.

  ! Attention: hydromet_list needs to be set up before this routine is called.
  !
  ! Bogus example
  ! References:
  !
  ! None
  !-----------------------------------------------------------------------


    use array_index, only: &
        hydromet_list  ! Names of the hydrometeor species

    use stats_variables, only: &
        irrm_bt,   & ! Variable(s)
        irrm_mc,   &
        irrm_hf,   &
        irrm_wvhf, &
        irrm_cl,   &
        irim_bt,   &
        irim_mc,   &
        irim_hf,   &
        irim_wvhf, &
        irim_cl,   &
        irgm_bt,   &
        irgm_mc,   &
        irgm_hf,   &
        irgm_wvhf, &
        irgm_cl,   &
        irsm_bt,   &
        irsm_mc,   &
        irsm_hf,   &
        irsm_wvhf, &
        irsm_cl

    use stats_variables, only: &
        iNrm_bt, & ! Variable(s)
        iNrm_mc, &
        iNrm_cl, &
        iNim_bt, &
        iNim_cl, &
        iNim_mc, &
        iNsm_bt, &
        iNsm_cl, &
        iNsm_mc, &
        iNgm_bt, &
        iNgm_cl, &
        iNgm_mc, &
        iNcm_bt, &
        iNcm_cl, &
        iNcm_mc

    use clubb_precision, only: &
        core_rknd

    use constants_clubb, only: &
        zero

    implicit none

    ! Input Variables
    integer, intent(in) :: ihm

    ! Input/Output Variables
    real( kind = core_rknd ), intent(inout) :: &
      max_velocity ! Maximum sedimentation velocity [m/s]

    integer, intent(inout) :: ixrm_hf, ixrm_wvhf, ixrm_cl, &
                              ixrm_bt, ixrm_mc

  !-----------------------------------------------------------------------

    !----- Begin Code -----

    ! Initializing max_velocity in order to avoid a compiler warning.
    ! Regardless of the case, it will be reset in the 'select case'
    ! statement immediately below.
    max_velocity = zero

    select case ( trim( hydromet_list(ihm) ) )
      case ( "rrm" )
        ixrm_bt   = irrm_bt
        ixrm_hf   = irrm_hf
        ixrm_wvhf = irrm_wvhf
        ixrm_cl   = irrm_cl
        ixrm_mc   = irrm_mc

        max_velocity = -9.1_core_rknd ! m/s

      case ( "rim" )
        ixrm_bt   = irim_bt
        ixrm_hf   = irim_hf
        ixrm_wvhf = irim_wvhf
        ixrm_cl   = irim_cl
        ixrm_mc   = irim_mc

        max_velocity = -1.2_core_rknd ! m/s

      case ( "rsm" )
        ixrm_bt   = irsm_bt
        ixrm_hf   = irsm_hf
        ixrm_wvhf = irsm_wvhf
        ixrm_cl   = irsm_cl
        ixrm_mc   = irsm_mc

        ! Morrison limit
!         max_velocity = -1.2_core_rknd ! m/s
        ! Made up limit.  The literature suggests that it is quite possible
        ! that snow flake might achieve a terminal velocity of 2 m/s, and this
        ! happens in the COAMPS microphysics -dschanen 29 Sept 2009
        max_velocity = -2.0_core_rknd ! m/s

      case ( "rgm" )
        ixrm_bt   = irgm_bt
        ixrm_hf   = irgm_hf
        ixrm_wvhf = irgm_wvhf
        ixrm_cl   = irgm_cl
        ixrm_mc   = irgm_mc

        max_velocity = -20._core_rknd ! m/s

      case ( "Nrm" )
        ixrm_bt   = iNrm_bt
        ixrm_hf   = 0
        ixrm_wvhf = 0
        ixrm_cl   = iNrm_cl
        ixrm_mc   = iNrm_mc

        max_velocity = -9.1_core_rknd ! m/s

      case ( "Nim" )
        ixrm_bt   = iNim_bt
        ixrm_hf   = 0
        ixrm_wvhf = 0
        ixrm_cl   = iNim_cl
        ixrm_mc   = iNim_mc

        max_velocity = -1.2_core_rknd ! m/s

      case ( "Nsm" )
        ixrm_bt   = iNsm_bt
        ixrm_hf   = 0
        ixrm_wvhf = 0
        ixrm_cl   = iNsm_cl
        ixrm_mc   = iNsm_mc

        ! Morrison limit
!         max_velocity = -1.2_core_rknd ! m/s
        ! Made up limit.  The literature suggests that it is quite possible
        ! that snow flake might achieve a terminal velocity of 2 m/s, and this
        ! happens in the COAMPS microphysics -dschanen 29 Sept 2009
        max_velocity = -2.0_core_rknd ! m/s

      case ( "Ngm" )
        ixrm_bt   = iNgm_bt
        ixrm_hf   = 0
        ixrm_wvhf = 0
        ixrm_cl   = iNgm_cl
        ixrm_mc   = iNgm_mc

        max_velocity = -20._core_rknd ! m/s

      case ( "Ncm" )
        ixrm_bt   = iNcm_bt
        ixrm_hf   = 0
        ixrm_wvhf = 0
        ixrm_cl   = iNcm_cl
        ixrm_mc   = iNcm_mc

        ! Use the rain water limit, since Morrison has no explicit limit on
        ! cloud water.  Presumably these numbers are never large.
        ! -dschanen 28 Sept 2009
        max_velocity = -9.1_core_rknd ! m/s

      case default
        ixrm_bt   = 0
        ixrm_hf   = 0
        ixrm_wvhf = 0
        ixrm_cl   = 0
        ixrm_mc   = 0

        max_velocity = -9.1_core_rknd ! m/s

    end select


    return

  end subroutine setup_stats_indices
  !-----------------------------------------------------------------------

end module fill_holes
